
import numpy as np
import matplotlib.pyplot as plt
from pylab import imshow,gray,show
from mpl_toolkits import mplot3d
import re
import matplotlib.cm as cm
from matplotlib.colors import Normalize
import matplotlib.pyplot as plotLine
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import LinearRegression
import statsmodels.api as sm
import statistics as stat
from sklearn.metrics import mean_squared_error
import itertools
import math
import matplotlib.transforms as mtransforms
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter
import os
from sklearn import preprocessing


def get_time_dipole(name, dcdfreq, flaginittime):
    data=np.loadtxt(f'{name}',dtype=float)
    
    #print(data)
    time=[]
    dipolez=[]
    anglex=[]
    angley=[]
    anglez=[]
    for row in data:
        if flaginittime==1:
            time.append((row[0]*dcdfreq*2e-6)+0.005)
            
        else: 
            time.append(row[0]*dcdfreq*2e-6)
            
        #time.append(row[0])
        #dipole.append(row[4])
        dipolez.append(row[3])
        anglex.append(math.degrees(math.acos(row[1]/row[4])))
        angley.append(math.degrees(math.acos(row[2]/row[4])))
        anglez.append(math.degrees(math.acos(row[3]/row[4])))
    mean_anglez=stat.mean(anglez)
    st_dev=stat.stdev(anglez)
    #print(name,"mean:", mean_anglez,"st_dev:",st_dev)
    return time, dipolez, anglex, angley, anglez, mean_anglez , st_dev

def get_average_anglez(data,init,finish,name):
    data_sel=[]
    #selet data 
    for i in range(init,finish+1):
        data_sel.append(data[i])
    mean_anglez=stat.mean(data_sel)
    st_dev=stat.stdev(data_sel)/(len(mean_anglez))**0.5
    print(name,"mean:", mean_anglez,"st_dev:",st_dev)
    return mean_anglez , st_dev
        
def get_average_data(Pz,theta,init,finish,name):
    Pz_sel=[]
    theta_sel=[]
    #selet data 
    for i in range(init,finish):
        Pz_sel.append(Pz[i])
        theta_sel.append(theta[i])
        
    mean_Pz=stat.mean(Pz_sel)
    st_dev_Pz=stat.stdev(Pz_sel)
    mean_theta=stat.mean(theta_sel)
    st_dev_theta=stat.stdev(theta_sel)
    print(name)
    print("mean_Pz:",mean_Pz,"st_dev_Pz:",st_dev_Pz)
    print("mean_theta:",mean_theta,"st_dev_theta:",st_dev_theta)
    num_bins=20
#     plt.hist(theta, num_bins)
#     plt.show()
#     plt.hist(Pz, num_bins)
#     plt.show()
    return mean_Pz, st_dev_Pz,mean_theta,st_dev_theta
        
    
def write_point(x, y, ax,hpos, vpos, xpos, ypos):
    x=round(x, 3)
    y=round(y, 3)
    trans_offset = mtransforms.offset_copy(ax.transData, fig=ax, x=xpos, y=ypos, units='dots')

    ax.text(x, y, f' {x}ns', verticalalignment=f'{vpos}', horizontalalignment=f'{hpos}',fontsize=18,zorder=2,transform=trans_offset) #put fontsize 11 por DW-0.2V/nm and another case 18


    ax.scatter(x, y, color='black', marker='o',zorder=4)

def get_time_energy(name, dcdfreq):
    data=np.loadtxt(f'{name}',dtype=float)
    #print(data)
    time=[]
    energy=[]
    for row in data:
        #time.append(row[0]*dcdfreq*2e-6)
        time.append((row[0]*dcdfreq*2e-6)+0.005)
        #time.append(row[0])
        energy.append(row[1])
    return time, energy

def compute_max_normalzD8C(field,name_sim):
    maxsPw=[]
    maxsPm=[]
    sim=[f'{location}{name_sim}_ZDMW.dat',f'{location}{name_sim}.dat']
    
    for i in field:
         sim.append (f'{location}{name_sim}0{i}E_ZDMW.dat')
         sim.append(f'{location}{name_sim}0{i}E.dat')
         
    for enum,n in enumerate(sim):
        if enum%2==0:
            maxsPw.append(max(get_time_dipole(n,2500,0)[1]))
        else:
            maxsPm.append(max(get_time_dipole(n,2500,0)[1]))
            
    Pmax=[max(maxsPw),max(maxsPm)]
    print("Pz max:", Pmax)
    return(Pmax)


location='/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/'
Fntsize=18

#plot settings 


fig = plt.figure()
fig.set_figheight(14)
fig.set_figwidth(20)
 
ax = plt.subplot2grid(shape=(2, 2), loc=(0, 0), colspan=2)
ax1 = plt.subplot2grid(shape=(2, 2), loc=(1, 0), colspan=1)
ax2 = plt.subplot2grid(shape=(2, 2), loc=(1, 1), colspan=1)


colors = cm.rainbow(np.linspace(0, 1, 7))
colors = np.arange(0, 1, 5)
colors=cm.viridis(np.linspace(0, 1, 7))
#files Z-DMw
DMPC=f'{location}DMPC.dat' #0V/nm
DMPCE2=f'{location}2DMPCE.dat' #0.1V/nm 
DMPCE0=f'{location}0DMPCE.dat' #0.2V/nm
#DMPCE0c=f'{location}0DMPCE.dat' #0.2V/nm
DMPCE4=f'{location}4DMPCE.dat' #0.3V/nm
DMPCE5=f'{location}5DMPCE.dat' #0.4V/nm
DMPCE1=f'{location}1DMPCE.dat' #0.5V/nm


#Z-DMw
ax.plot(get_time_dipole(DMPC,50000,1)[0],get_time_dipole(DMPC,50000,1)[1],label='DW-0V/nm',color=colors[0])
ax.plot(get_time_dipole(DMPCE2,2500,1)[0],get_time_dipole(DMPCE2,2500,1)[1],label='DW-0.1V/nm',color=colors[1])
#ax.plot(get_time_dipole(DMPCE3,2500,1)[0],get_time_dipole(DMPCE3,2500,1)[1],label='DW-0.16V/nm',color=colors[2])
ax.plot(get_time_dipole(DMPCE0,2500,1)[0],get_time_dipole(DMPCE0,2500,1)[1],label='DW-0.2V/nm',color=colors[2])
#ax.plot(get_time_dipole(DMPCE0c,2500,1)[0],get_time_dipole(DMPCE0c,2500,1)[1],label='DWc-0.2V/nm',color='black')
ax.plot(get_time_dipole(DMPCE4,2500,1)[0],get_time_dipole(DMPCE4,2500,1)[1],label='DW-0.3V/nm',color=colors[3])
ax.plot(get_time_dipole(DMPCE5,2500,1)[0],get_time_dipole(DMPCE5,2500,1)[1],label='DW-0.4V/nm',color=colors[4])
ax.plot(get_time_dipole(DMPCE1,2500,1)[0],get_time_dipole(DMPCE1,2500,1)[1],label='DW-0.5V/nm',color=colors[5])


#angle of the dipole moment with respect to the z plane.
ax1.plot(get_time_dipole(DMPC,50000,0)[0],get_time_dipole(DMPC,50000,1)[4],label='DW-0V/nm',color=colors[0])
ax1.plot(get_time_dipole(DMPCE2,2500,1)[0],get_time_dipole(DMPCE2,2500,1)[4],label='DW-0.1V/nm',color=colors[1])
#ax1.plot(get_time_dipole(DMPCE3,2500,1)[0],get_time_dipole(DMPCE3,2500,1)[4],label='DW-0.16V/nm',color=colors[2])
ax1.plot(get_time_dipole(DMPCE0,2500,1)[0],get_time_dipole(DMPCE0,2500,1)[4],label='DW-0.2V/nm',color=colors[2])
#x1.plot(get_time_dipole(DMPCE0c,2500,1)[0],get_time_dipole(DMPCE0c,2500,1)[4],label='DWc-0.2V/nm',color='black')
ax1.plot(get_time_dipole(DMPCE4,2500,1)[0],get_time_dipole(DMPCE4,2500,1)[4],label='DW-0.3V/nm',color=colors[3])
ax1.plot(get_time_dipole(DMPCE5,2500,1)[0],get_time_dipole(DMPCE5,2500,1)[4],label='DW-0.4V/nm',color=colors[4])
ax1.plot(get_time_dipole(DMPCE1,2500,1)[0],get_time_dipole(DMPCE1,2500,1)[4],label='DW-0.5V/nm',color=colors[5])


#plot names and legend
ax.set_xlabel('Time(ns)',fontsize=Fntsize)
ax.set_ylabel('$P_{zDMPC}$ (D)',fontsize=Fntsize)
ax1.set_ylabel('$\Theta_{zDMPC}$ (Degree)',fontsize=Fntsize)
ax1.set_xlabel('Time(ns)',fontsize=Fntsize)
ax.set_title('A)',fontsize=Fntsize)
ax1.set_title('B)',fontsize=Fntsize)
ax1.legend(loc="upper right",fontsize=Fntsize,ncol=2)
#ax.legend(loc="center right")
plt.tight_layout(pad=2)
ax.set_xlim(0,4)
ax1.set_xlim(0,4)

tp=[0,191,36,25,5]
#Put points 

#0.2V/nm
#write_point((get_time_dipole(DMPCE0,2500,1)[0])[145],(get_time_dipole(DMPCE0,2500,1)[1])[125],ax,'right')
write_point((get_time_dipole(DMPCE0,2500,1)[0])[191],(get_time_dipole(DMPCE0,2500,1)[1])[191],ax,'left','center',-15,13)

#0.3V/nm
write_point((get_time_dipole(DMPCE4,2500,1)[0])[36],(get_time_dipole(DMPCE4,2500,1)[1])[36],ax,'left','center',5,0)

#0.4V/nm
write_point((get_time_dipole(DMPCE5,2500,1)[0])[25 ],(get_time_dipole(DMPCE5,2500,1)[1])[25],ax,'left','center',-4,-10)

#0.5 v/nm
write_point((get_time_dipole(DMPCE1,2500,1)[0])[5],(get_time_dipole(DMPCE1,2500,1)[1])[5],ax,'left', 'center',-8, 24)



tp=[(get_time_dipole(DMPCE0,2500,1)[0])[191],(get_time_dipole(DMPCE4,2500,1)[0])[36],(get_time_dipole(DMPCE5,2500,1)[0])[25],(get_time_dipole(DMPCE1,2500,1)[0])[5]]
E=[0.2,0.3,0.4,0.5]

ax2.scatter(E,tp)
ax2.plot(E,tp, color='black')
ax2.set_title('C)',fontsize=Fntsize)
ax2.set_ylabel('$t_{p} (ns)$',fontsize=Fntsize)
ax2.set_xlabel('$E$ $(V/nm)$',fontsize=Fntsize)
plt.tight_layout()

plt.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/Z-DMW_DW.png')

print((tp[0]-tp[1]))
print((tp[1]-tp[2]))

0.7749999999999999
0.05499999999999999


AvgDM_DMPC_0={}
AvgDM_Water_0={}
AvgDM_DMPC_01={}
AvgDM_Water_01={}


# fig, axs = plt.subplots(1, 1,figsize=(16, 8))

DMPCE0=f'{location}0DMPCE_ZDMW.dat' #0.2V/nm
DMPCE0o=f'{location}0DMPCE_ZDMW_origin.dat' #0.2V/nm
DMPCE0e=f'{location}0DMPCE_ZDMW_e.dat' #0.2V/nm
# ax=axs




fig, axes = plt.subplots(3, 4, figsize=(10, 10))
fig.set_figheight(14)
fig.set_figwidth(15)
axes = axes.flatten()
axes[0] = plt.subplot2grid(shape=(3, 4), loc=(0, 0), colspan=3)

#Compute derivate of the dipole moment
xvalues=get_time_dipole(DMPCE0,2500,1)[0]
yvalues=get_time_dipole(DMPCE0,2500,1)[1]

#xvalues = savgol_filter(xvalues, 9, 2)
yvalues = savgol_filter(yvalues, 41, 4)
dZ_DMW=np.gradient(yvalues,xvalues)
axes[0].plot(xvalues,dZ_DMW,label=r'$\frac{dP_{zw}}{dt}$',color='gray')
axes[0].plot(get_time_dipole(DMPCE0,2500,1)[0],get_time_dipole(DMPCE0,2500,1)[1],label=r'$P_{zw}$',color=colors[2])

axes[0].set_xlabel('Time(ns)',fontsize=Fntsize)
axes[0].set_ylabel(r'$P_{zw} (D) \ , \  \frac{dP_{zw}}{dt} (D/ns)$',fontsize=Fntsize)
ax=axes[0]
ax.set_xlim(0,4)
#ax.set_ylim(-50000,50000)
ax.legend(loc="upper right",fontsize=Fntsize,ncol=2)
times=[0,122,191,241,268,304,332,441,799]
a=3
for t in times:
    img=plt.imread(f'/media/thalia/New Volume/SIMULACIONES TDG/images TDG/DMPCE/0DMPCE/0DMPCE_{t}.png')
    axes[a].imshow(img)
    #axes[a].axis('off')  # Turn off axis labels
    tns=round((t*2500*2e-6)+0.005, 3)
    axes[a].set_title(f'{tns}ns',fontsize=Fntsize)
    axes[a].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)   
    axes[a].xaxis.set_visible(False)
    #axes[4].yaxis.set_visible(False)
    axes[3].set_ylabel('E',fontsize=Fntsize)
    axes[a].tick_params(left=False, labelleft=False)
    plt.setp(axes[a].spines.values(), visible=False)
    a+=1
start_point = (0, 556)
end_point = (0, 100)
axes[3].arrow(start_point[0], start_point[1], 
          end_point[0] - start_point[0], end_point[1] - start_point[1],
          color='black', width=3, head_width=10, head_length=10)
axes[3].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)   
axes[3].xaxis.set_visible(False)
#axes[4].yaxis.set_visible(False)
axes[3].set_ylabel('$E_{z}=0.2V/nm$',fontsize=20)
axes[3].tick_params(left=False, labelleft=False)
plt.setp(axes[3].spines.values(), visible=False)


# ax.plot(get_time_dipole(DMPCE0o,2500,1)[0],get_time_dipole(DMPCE0o,2500,1)[1],label='DW-0.2V/nm_origin',color=colors[5])
# ax.plot(get_time_dipole(DMPCE0e,2500,1)[0],get_time_dipole(DMPCE0e,2500,1)[1],label='DW-0.2V/nm_e',color=colors[6])


#Put points 

#0.2V/nm
write_point((get_time_dipole(DMPCE0,2500,1)[0])[0],(get_time_dipole(DMPCE0,2500,1)[1])[0],axes[0],'right','center',33,-7)
write_point((get_time_dipole(DMPCE0,2500,1)[0])[122],(get_time_dipole(DMPCE0,2500,1)[1])[122],ax,'right','center',10,-13)
write_point((get_time_dipole(DMPCE0,2500,1)[0])[191],(get_time_dipole(DMPCE0,2500,1)[1])[191],ax,'right','center',10,-13)
write_point((get_time_dipole(DMPCE0,2500,1)[0])[241],(get_time_dipole(DMPCE0,2500,1)[1])[241],ax,'right','center',10,-10)
write_point((get_time_dipole(DMPCE0,2500,1)[0])[268],(get_time_dipole(DMPCE0,2500,1)[1])[268],ax,'left','center',3,-4)
write_point((get_time_dipole(DMPCE0,2500,1)[0])[304],(get_time_dipole(DMPCE0,2500,1)[1])[304],ax,'left','center',5,0)
write_point((get_time_dipole(DMPCE0,2500,1)[0])[332],(get_time_dipole(DMPCE0,2500,1)[1])[332],ax,'left','center',5,0)
write_point((get_time_dipole(DMPCE0,2500,1)[0])[441],(get_time_dipole(DMPCE0,2500,1)[1])[441],ax,'left','center',1,8)
write_point((get_time_dipole(DMPCE0,2500,1)[0])[799],(get_time_dipole(DMPCE0,2500,1)[1])[799],ax,'left','center',-30,-10)
#Adjust layout to prevent overlapping
plt.tight_layout()
plt.show()
show()
fig.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/Z-DMW_DW_0.2.png',dpi=300)


from PIL import Image
def conv_frame_ns(frames,dcdfreq):
    timesns=[]
    for f in frames:
      timesns.append(round((f*dcdfreq*2e-6)+0.005,2))
    return(timesns)


loctp='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/tps/'
loctp20='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/tps20/'
loctpend='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/the end/'
loctpendtop='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/theEndtop/'
locations=[loctp,loctp20,loctpend,loctpendtop]
#psimage.save('myImage.png')
# fig, axs =pl


fig, axes = plt.subplots(4, 4, figsize=(10, 10))
fig.set_figheight(20)
fig.set_figwidth(20)
row=0
tpDMPCE0=[191,46,35,5]
simulations=[0,4,5,1]
#simulations=[0]
timetpns=conv_frame_ns(tpDMPCE0,2500)
timetp20ns=conv_frame_ns([fr+20 for fr in tpDMPCE0],2500)
timetpendns=conv_frame_ns([799,799,231,253],2500)
alltimens=[timetpns,timetp20ns,timetpendns,timetpendns]
fieldname=['0.2V/nm','0.3V/nm','0.4V/nm','0.5V/nm']
z='z'
start_point = (-1, 356)
end_point = (-1, 100)
name='DMPCE'
for row,loc in enumerate(locations):
    timesns=alltimens[row]
    for col,idx in enumerate(simulations):

        psimage=Image.open(f'{loc}{idx}{name}')
        axes[row,col].imshow(psimage)
        axes[row,col].set_title(f'{timesns[col]}ns',fontsize=Fntsize+5)
        axes[row,col].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)   
        axes[row,col].xaxis.set_visible(False)
        axes[row,col].tick_params(left=False, labelleft=False)
        plt.setp(axes[row,col].spines.values(), visible=False)
        #axes[row,col].yaxis.set_visible(False)
        #plot field vector
        if row==0:
            axes[row,col].arrow(start_point[0], start_point[1], 
                      end_point[0] - start_point[0], end_point[1] - start_point[1],
                      color='black', width=6, head_width=15, head_length=10)
        
            axes[row,col].set_ylabel(f'$E_{z}$={fieldname[col]}',fontsize=20)
fig.suptitle('DW',fontsize=30)
# #plt.subplots_adjust(left=0.1,
#                     bottom=0.1, 
#                     right=0.9, 
#                     top=0.9, 
#                     wspace=0.1, 
#                     hspace=0)       
plt.tight_layout(pad=0.5)

fig.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/timetp_DW.png',dpi=300)


from PIL import Image
def conv_frame_ns(frames,dcdfreq):
    timesns=[]
    for f in frames:
      timesns.append(round((f*dcdfreq*2e-6)+0.005,2))
    return(timesns)

loctin='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/DMin/'
loctp='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/DMtps/'
loctp20='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/DMtps20/'
loctpend='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/DMend/'
loctpendtop='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/theEndtop/'
locations=[loctin, loctp,loctp20,loctpend]
#psimage.save('myImage.png')
# fig, axs =pl


fig, axes = plt.subplots(2, 4, figsize=(10, 10))
fig.set_figheight(10)
fig.set_figwidth(15)
row=0
tp=[191,5]
simulations=[0,1]
#simulations=[0]
timetpns=conv_frame_ns(tp,2500)
timetp20ns=conv_frame_ns([fr+20 for fr in tp],2500)
timetpendns=conv_frame_ns([799,253],2500)
alltimens=[[0.005]*8,timetpns,timetp20ns,timetpendns,timetpendns]
fieldname=['0.2V/nm','0.5V/nm']
z='z'
start_point = (-1, 356)
end_point = (-1, 100)
name='DMPCE'
for row,loc in enumerate(locations):
    timesns=alltimens[row]
    for col,idx in enumerate(simulations):

        psimage=Image.open(f'{loc}{idx}{name}')
        axes[col,row].imshow(psimage)
        axes[col,row].set_title(f'{timesns[col]}ns',fontsize=Fntsize+5)
        axes[col,row].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)   
        axes[col,row].xaxis.set_visible(False)
        axes[col,row].tick_params(left=False, labelleft=False)
        plt.setp(axes[col,row].spines.values(), visible=False)
        if row==0:
            axes[col,row].arrow(start_point[0], start_point[1], 
                      end_point[0] - start_point[0], end_point[1] - start_point[1],
                      color='black', width=6, head_width=20, head_length=10)
        
            axes[col,row].set_ylabel(f'$E_{z}$={fieldname[col]}',fontsize=20)
fig.suptitle('DW',fontsize=30)
     
plt.tight_layout(pad=0.5)

fig.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/Pdmpcarrowtimetp.png',dpi=300)


from PIL import Image
def conv_frame_ns(frames,dcdfreq):
    timesns=[]
    for f in frames:
      timesns.append(round((f*dcdfreq*2e-6)+0.005,3))
    print(timesns)
    return(timesns)

loc='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/DMlipid/'
names=['05lp0','05lp40','05lp253']
#psimage.save('myImage.png')
# fig, axs =pl


fig, axes = plt.subplots(1, 4, figsize=(10, 10))
fig.set_figheight(7)
fig.set_figwidth(18)
axes.flatten()
row=0
tp=[0,40, 253]
sim='/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/lipid24_05.dat'
data=get_time_dipole(sim,2500,1)
t=data[0]
theta=data[4]
axes[0].plot(t,theta)
axes[0].set_xlabel('Time(ns)',fontsize=Fntsize+4)
axes[0].set_ylabel('$\Theta_{z}$ (Degree)',fontsize=Fntsize+4)
timetpns=conv_frame_ns(tp,2500)
print(timetpns)
alltimens=timetpns
fieldname=['0.2V/nm','0.5V/nm']
z='z'
start_point = (-1, 356)
end_point = (-1, 100)

for row, name in enumerate(names):
        row +=1
        psimage=Image.open(f'{loc}{name}')
        axes[row].imshow(psimage)
        axes[row].set_title(f'{alltimens[row-1]}ns',fontsize=Fntsize+5)
        axes[row].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)   
        axes[row].xaxis.set_visible(False)
        axes[row].tick_params(left=False, labelleft=False)
        plt.setp(axes[row].spines.values(), visible=False)

        if row==1:
            axes[row].arrow(start_point[0], start_point[1], 
                      end_point[0] - start_point[0], end_point[1] - start_point[1],
                      color='black', width=6, head_width=20, head_length=10)
        
            axes[row].set_ylabel(f'$E_{z}$={fieldname[row]}',fontsize=20)
fig.suptitle('DW: 0.5V/nm',fontsize=30)
  
plt.tight_layout(pad=0.5)

fig.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/lipidDM.png',dpi=300)

[0.005, 0.205, 1.27]
[0.005, 0.205, 1.27]


from PIL import Image
def conv_frame_ns(frames,dcdfreq):
    timesns=[]
    for f in frames:
      timesns.append(round((f*dcdfreq*2e-6)+0.005,3))
    print(timesns)
    return(timesns)

loc='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/Effectpepion/'
names=['DW05','D8CV02E','1DI02E']
#psimage.save('myImage.png')
# fig, axs =pl


fig, axes = plt.subplots(3, 4, figsize=(10, 10))
fig.set_figheight(18)
fig.set_figwidth(18)
#axes.flatten()
row=0
tp=[0,40, 253]

timetpns=conv_frame_ns(tp,2500)
print(timetpns)
alltimens=timetpns
fieldname=['0.5V/nm','0.2V/nm','0.2V/nm','0.2V/nm']
simname=['DW','D8CV', '1DI']
z='z'
start_point = (-1, 356)
end_point = (-1, 100)

times=[0,99,230,230 ]
alltimens=conv_frame_ns(times,2500)

for row, name in enumerate(names):
        for col, frame in enumerate([0,99,230,230]): 
            
            if frame==230 and col==3:
                psimage=Image.open(f'{loc}{name}_{frame}top')
                print("hola")
            else:

                psimage=Image.open(f'{loc}{name}_{frame}')
                
            axes[row,col].imshow(psimage)
            if col==0 and row==0:
                axes[row,col].set_title(f'{simname[row]}:   {alltimens[col]}ns',fontsize=Fntsize+5)
            elif col==0:
                axes[row,col].set_title(f'{simname[row]}', fontsize=Fntsize+5)
            elif row==0:
                axes[row,col].set_title(f'{alltimens[col]}ns',fontsize=Fntsize+5)
            axes[row,col].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)   
            axes[row,col].xaxis.set_visible(False)
            axes[row,col].tick_params(left=False, labelleft=False)
            plt.setp(axes[row,col].spines.values(), visible=False)
            
            if col==0:
                axes[row,col].arrow(start_point[0], start_point[1], 
                          end_point[0] - start_point[0], end_point[1] - start_point[1],
                          color='black', width=6, head_width=20, head_length=10)

                axes[row,col].set_ylabel(f'$E_{z}$={fieldname[row]}',fontsize=20)
     
plt.tight_layout(pad=0.5)

fig.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/effecpepion.png',dpi=300)

[0.005, 0.205, 1.27]
[0.005, 0.205, 1.27]
[0.005, 0.5, 1.155, 1.155]
hola
hola
hola


##########
## D8CV ##
##########  
plt.rc('font', size=18)
fig, ax = plt.subplots(2, 2,figsize=(25, 15))

# fig = plt.figure()
# fig.set_figheight(14)
# fig.set_figwidth(20)
step=101
odr=9
ypt=0.9
Fntsize=30
lw=3
location='/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/'
#simulation parameter
dif_con_sim=[4,2,5,3] #[4,1,6,3]
name_sim='0D8CV'
#Pmax=compute_max_normalzD8C(range(1,6),name_sim)

#Smoothing parameter


#plot parameters
transparency=0.3 #transparency of the data without simuation
#lw=3 #line with of the plot
labels=['0.0V/nm']
titles=["$E_{z}$", "-$E_{z}$"]
colors=['black']
lines=['solid']
#ax3 = ax.twinx()

#plot D8CV 0.0V/nm
for i,sim in enumerate([f'{location}{name_sim}_ZDMW.dat',f'{location}{name_sim}.dat']):
        s=0
        lw=2
        x=np.array(get_time_dipole(sim,2500,1)[0])
        y=np.array(get_time_dipole(sim,2500,1)[1])
        theta=get_time_dipole(sim,2500,1)[4]


        #suavizado
        theta_smt=savgol_filter(theta,step,odr)
        y_smt=savgol_filter(y,step,odr)
    #             ynorm_smt=y_smt/Pmax[i]
    #             ynorm=y/Pmax[i]
        ynorm_smt=y_smt
        ynorm=y
        
        if i==0 : 
            ax[0,0].plot(x,ynorm_smt,color=colors[s],linestyle=lines[s], label=f'{labels[s]}',linewidth=lw)
            ax[0,0].plot(x,ynorm,color=colors[s],linestyle=lines[s], alpha=0.3)
            ax[0,1].plot(x,theta_smt,color=colors[s],label=f'{labels[s]}',linestyle=lines[s],linewidth=lw)

                
                #axes[2,2].legend(bbox_to_anchor=(2.5,-0.2),fontsize=Fntsize+5,ncol=9) 
           #plot Pz and theta of DMPC
        else:
            ax[1,0].plot(x,ynorm_smt,color=colors[s],linestyle=lines[s],linewidth=lw,label=f'{labels[s]}')
            ax[1,0].plot(x,ynorm,color=colors[s],linestyle=lines[s], alpha=0.3)
            ax[1,1].plot(x,theta_smt,color=colors[s],linewidth=lw,linestyle=lines[s],label=f'{labels[s]}')

labels=[['0.1V/nm', '0.2V/nm'],['-0.1V/nm', '-0.2V/nm']]
colors=cm.rainbow(np.linspace(0, 1, 4))
lw=3
fields=[1,2]
tp=[99,68]
for row, namesimsufix in enumerate(['E','Er']):
    for s , E in enumerate(fields):
        
        #print("s:",s)
        if E==1:
            simulations=[f'{location}{name_sim}0{E}{namesimsufix}_ZDMW.dat',f'{location}{name_sim}0{E}{namesimsufix}_DMPC.dat']

        else:
            simulations=[f'{location}{name_sim}0{E}{namesimsufix}_ZDMW.dat',f'{location}{name_sim}0{E}{namesimsufix}.dat']
        print(simulations)
        
        lines=['solid','dotted']
        a=0

        for i in range(0,len(simulations)): 
            sim=simulations[i] #0.1V/nm
            x=np.array(get_time_dipole(sim,2500,1)[0])
            y=np.array(get_time_dipole(sim,2500,1)[1])
            theta=get_time_dipole(sim,2500,1)[4]


            #suavizado
            theta_smt=savgol_filter(theta,step,odr)
            y_smt=savgol_filter(y,step,odr)
        #             ynorm_smt=y_smt/Pmax[i]
        #             ynorm=y/Pmax[i]
            ynorm_smt=y_smt
            ynorm=y

            #plot Pz and theta of water
            if i==0 : 
                ax[0,0].plot(x,ynorm_smt,color=colors[s],linestyle=lines[row], label=f'{labels[row][s]}',linewidth=lw)
                ax[0,0].plot(x,ynorm,color=colors[s],linestyle=lines[row], alpha=0.3)
                ax[0,1].plot(x,theta_smt,color=colors[s],label=f'{labels[row][s]}',linestyle=lines[row],linewidth=lw)
                
                
                #axes[2,2].legend(bbox_to_anchor=(2.5,-0.2),fontsize=Fntsize+5,ncol=9) 
           #plot Pz and theta of DMPC
            else:
                ax[1,0].plot(x,ynorm_smt,color=colors[s],linestyle=lines[row],linewidth=lw,label=f'{labels[row][s]}')
                ax[1,0].plot(x,ynorm,color=colors[s],linestyle=lines[row], alpha=0.3)
                ax[1,1].plot(x,theta_smt,color=colors[s],linewidth=lw,linestyle=lines[row],label=f'{labels[row][s]}')

                
            if s==1 and i==1:
                ax[0,0].axvline(x=x[tp[row]], color=colors[s], linestyle=lines[row])
                ax[0,1].axvline(x=x[tp[row]], color=colors[s], linestyle=lines[row])
                ax[1,0].axvline(x=x[tp[row]], color=colors[s], linestyle=lines[row])
                ax[1,1].axvline(x=x[tp[row]], color=colors[s], linestyle=lines[row])
#         ax[row,s].set_xlim(0,4) 

    ax[row,1].set_ylim(0,180)
    

#plot dipole:

from PIL import Image
def conv_frame_ns(frames,dcdfreq):
    timesns=[]
    for f in frames:
      timesns.append(round((f*dcdfreq*2e-6)+0.005,2))
    return(timesns)




loc='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/dipoleMomentumDC/'
location='/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/'


           
ax[0,0].set_ylabel('$P_{zw} (D)$',fontsize=Fntsize+5)
 
ax[1,0].set_ylabel('$P_{zDMPC} (D)$',fontsize=Fntsize+5)
#ax[2,0].set_ylabel('$P_{zall-pep} (D)$',fontsize=Fntsize+5)
ax[0,1].set_ylabel('$\Theta_{zw}$ (Degree)',fontsize=Fntsize+5)  
ax[1,1].set_ylabel('$\Theta_{zDMPC}$ (Degree)',fontsize=Fntsize+5)
#ax[2,1].set_ylabel('$\Theta_{zall-pep}$ (Degree)',fontsize=Fntsize+5)
ax[1,0].set_xlim(0,6) 
ax[1,1].set_xlim(0,6) 
ax[1,0].set_xlabel('Time(ns)',fontsize=Fntsize+5)
ax[1,1].set_xlabel('Time(ns)',fontsize=Fntsize+5)
ax[1,0].legend(bbox_to_anchor=(2.1, -0.19),fontsize=Fntsize,ncol=6) 
#plt.tight_layout(pad=3)
plt.show()
plt.show()
fig.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/P_W_dmcp_p_D8CV.png',dpi=300) 



chainpep={"DCV": ['A'], "D3CV": ['A','B','C'], "D8CV": ['A','B','C','D','E','F','G','H'],"DCH": ['A'], "D3CH": ['A','B','C']}
fig, axes = plt.subplots(4, 4, figsize=(10, 10))
fig.set_figheight(35)
fig.set_figwidth(40)
names_snap=['0D8CV01E','0D8CV01Er','0D8CV02E','0D8CV02Er']
titles=['D8CV: 0.1V/nm','D8CV: -0.1V/nm','D8CV: 0.2V/nm','D8CV: -0.2V/nm' ] 
tppnD8CV=[99,68]
colors=['yellow','gray', 'orange','red', 'magenta', 'pink', 'green','indigo', 'black']




for col,name in enumerate(names_snap):
    
    psimage=Image.open(f'{loc}{name}')
    axes[0,col].imshow(psimage)
    axes[0,col].set_title(f'{titles[col]}',fontsize=Fntsize+5)
    axes[0,col].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)   
    axes[0,col].xaxis.set_visible(False)
    axes[0,col].tick_params(left=False, labelleft=False)
    plt.setp(axes[0,col].spines.values(), visible=False)
    #axes[row,col].yaxis.set_visible(False)
    
    sufixs=chainpep['D8CV']
    for c,sufix in enumerate(sufixs): 
        
        sim=f'{location}{name}_p{sufix}.dat'
        data=get_time_dipole(sim,2500,1)
        t=data[0]
        P=data[1]
        theta=data[4]
        
        
        axes[2,col].plot(t,P,label=f'Pep_{sufix}',color=colors[c])
        axes[3,col].plot(t,theta,label=f'Pep{sufix}',color=colors[c])
        
names_snap=['0D8CV01Etop','0D8CV01Ertop','0D8CV02Etop','0D8CV02Ertop']        
for col,name in enumerate(names_snap):
    
    psimage=Image.open(f'{loc}{name}')
    axes[1,col].imshow(psimage)
    #axes[1,col].set_title(f'{titles[col]}',fontsize=Fntsize+5)
    axes[1,col].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)   
    axes[1,col].xaxis.set_visible(False)
    axes[1,col].tick_params(left=False, labelleft=False)
    plt.setp(axes[1,col].spines.values(), visible=False)

axes[2,0].set_ylabel('$P_{zp} (D)$',fontsize=Fntsize+5)
axes[3,0].set_ylabel('$\Theta_{zp}$ (Degree)',fontsize=Fntsize+5)
for i in range(0,4):
    axes[2,i].set_xlabel('Time(ns)',fontsize=Fntsize+5)


axes[3,2].legend(bbox_to_anchor=(2.4,-0.2),fontsize=Fntsize+5,ncol=9) 


plt.show()
fig.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/P_pnD8CV.png',dpi=300)

['/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/0D8CV01E_ZDMW.dat', '/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/0D8CV01E_DMPC.dat']
['/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/0D8CV02E_ZDMW.dat', '/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/0D8CV02E.dat']
['/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/0D8CV01Er_ZDMW.dat', '/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/0D8CV01Er_DMPC.dat']
['/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/0D8CV02Er_ZDMW.dat', '/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/0D8CV02Er.dat']


from PIL import Image
def conv_frame_ns(frames,dcdfreq):
    timesns=[]
    for f in frames:
      timesns.append(round((f*dcdfreq*2e-6)+0.005,2))
    return(timesns)

plt.rc('font', size=15) 

loc='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/dipoleMomentumDC/'
location='/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/'
names_snap=['DCV','D3CV','D8CV','DCH','D3CH']
chainpep={"DCV": ['A'], "D3CV": ['A','B','C'], "D8CV": ['A','B','C','D','E','F','G','H'],"DCH": ['A'], "D3CH": ['A','B','C']}
colors=['yellow','gray', 'orange','red', 'magenta', 'pink', 'green','indigo' ]

fig, axes = plt.subplots(3, 5, figsize=(10, 10))
fig.set_figheight(20)
fig.set_figwidth(35)

for col,name in enumerate(names_snap):
    
    psimage=Image.open(f'{loc}{name}')
    axes[0,col].imshow(psimage)
    axes[0,col].set_title(f'{name}',fontsize=Fntsize+5)
    axes[0,col].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)   
    axes[0,col].xaxis.set_visible(False)
    axes[0,col].tick_params(left=False, labelleft=False)
    plt.setp(axes[0,col].spines.values(), visible=False)
    #axes[row,col].yaxis.set_visible(False)
    
    sufixs=chainpep[name]
    for c,sufix in enumerate(sufixs): 
        
        sim=f'{location}{name}_p{sufix}.dat'
        data=get_time_dipole(sim,2500,1)
        t=data[0]
        P=data[1]
        theta=data[4]
        axes[1,col].plot(t,P,label=f'Pep_{sufix}',color=colors[c])
        axes[2,col].plot(t,theta,label=f'Pep_{sufix}',color=colors[c])
axes[1,0].set_ylabel('$P_{zp} (D)$',fontsize=Fntsize+5)
axes[2,0].set_ylabel('$\Theta_{zp}$ (Degree)',fontsize=Fntsize+5)
axes[2,2].set_xlabel('Time(ns)',fontsize=Fntsize+5)
axes[2,2].legend(bbox_to_anchor=(4,-0.2),fontsize=Fntsize+5,ncol=8) 
plt.show()
    
# plt.tight_layout(pad=0.5)

fig.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/pepDipoleDC.png',dpi=300)

Functions¶

Plot the behavior of the z-component of the membrane dipole moment and visualize the snapshots of the simulations.¶

Snapshot at $t_p$ and at the end of the simulation of the system consisting of water and the bilayer.¶

Compare P_DMPC and Pw of D8CVr¶

DIPOle momentum of peptides¶

Functions¶

Plot the behavior of the z-component of the membrane dipole moment and visualize the snapshots of the simulations.¶

Snapshot at tptpt_p and at the end of the simulation of the system consisting of water and the bilayer.¶

Compare P_DMPC and Pw of D8CVr¶

DIPOle momentum of peptides¶

Snapshot at $t_p$ and at the end of the simulation of the system consisting of water and the bilayer.¶